!             Program block renorm-vlowk.f90
!
!             Authors:  Morten Hjorth-Jensen
!             ADDRESS:  Dept. Physics, University Oslo, N-0316 OSLO
!             E-MAIL:   morten.hjorth-jensen@fys.uio.no
!             LANGUAGE: F90/F95  
!             LAST UPGRADE : March 2007
!
!                        
!
!                    Begin setup of similarity transformation
!
SUBROUTINE  setup_vkk
  USE partial_waves
  USE constants
  USE relcm_gmatrix
  USE single_particle_orbits
  USE wave_functions
  USE configurations
  IMPLICIT NONE
  INTEGER :: loop, number_orbits, max_conf,ncoup

  !     reserve space in memory for mesh point and h.o. wave function  arrays
  ALLOCATE ( ra (n_rel), wra (n_rel));   ALLOCATE ( krel (n_rel), wkrel (n_rel))
!     set up mesh points for relative coordinates
  CALL vlowk_mesh
  
  allocate( x_mesh(qint_points), wx_mesh(qint_points) )
  call q_mesh
  
! ALLOCATE( t_mesh (qint_points), wt_mesh (qint_points))
! call cheb_mesh

  ALLOCATE ( q_x(qint_points), q_w(qint_points) )

  
  !     find all partial waves with given jmin and jmax
  CALL setup_channels
  if (type_of_renormv =='vbare') then
     vkk_dim=n_k
  elseif (type_of_renormv =='vlowk') then
     vkk_dim=n_k1
  endif
  ALLOCATE(vfree(2*vkk_dim,2*vkk_dim,no_channels))

  DO loop = 1, no_channels
     ncoup=1
     IF ( orb_lrel_max(loop) == jang_rel(loop) ) ncoup = 1
     IF ( orb_lrel_max(loop) /= jang_rel(loop) ) ncoup = 2
     CALL vkk_channel(loop,ncoup)
  ENDDO

END SUBROUTINE setup_vkk

SUBROUTINE dealloc_vkk
  USE wave_functions
  USE relcm_gmatrix
  DEALLOCATE ( x_mesh,wx_mesh)
!  DEALLOCATE (t_mesh,wt_mesh)
  DEALLOCATE (vfree)
  DEALLOCATE ( q_x, q_w)
  DEALLOCATE ( ra, wra, krel, wkrel)  
END SUBROUTINE dealloc_vkk


! for test
subroutine v_ho(v_pot,na,la,ja,nb,lb,jb,nc,lc,jc,nd,ld,jd,j_lab,isospin_tz)
  use wave_functions
  implicit none
  REAL(DP),INTENT(INOUT) :: v_pot
  REAL(DP) :: hoa,hob,hoc,hod,v_t,pa,pb,pc,pd
  INTEGER, INTENT(IN) :: isospin_tz, j_lab, ja, jb, jc, jd, la, lb, lc, ld,na,nb,nc,nd
  INTEGER::i,j,k,l,iph
  v_pot=0.d0
  do i=1,n_lab
     do j=1,n_lab
        do k=1,n_lab
           do l=1,n_lab
              hoa=rnlr(i,la,na)
              hob=rnlr(j,lb,nb)
              hoc=rnlr(k,lc,nc)
              hod=rnlr(l,ld,nd)
              pa=k_mesh(i)
              pb=k_mesh(j)
              pc=k_mesh(k)
              pd=k_mesh(l)
              CALL v_k(v_t,la,ja,pa,lb,jb,pb,lc,jc,pc,ld,jd,pd,isospin_tz,j_lab)
              v_pot=v_pot+v_t*hoa*hob*hoc*hod
           enddo
        enddo
     enddo
  enddo
  v_pot=v_pot*iph(ABS(la+lb-lc-ld)/2)
end subroutine v_ho




!! Mev*fm^2, v_pot=<ka la ja kb lb jb|V|kc lc jc kd ld jd>, note that |n l> = \int dk R_{nl}(k) k |k l>, |vec(k)>=sum_{lm} |kl>/k * Y_{lm}^*
!! pa,pb,pc,pd: fm^(-1)
SUBROUTINE v_k(v_pot,la,ja,pa,lb,jb,pb,lc,jc,pc,   &
     ld,jd,pd,isospin_tz,j_lab)
  USE configurations
  USE constants
  USE single_particle_orbits
  USE partial_waves
  USE relcm_gmatrix
  USE wave_functions
  USE ang_mom_functions
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: isospin_tz, j_lab, ja, jb, jc, jd, la, lb, lc, ld
  INTEGER :: llmin, llmax, lambda_min, lambda_max, ncoup, channel,dim
  INTEGER :: spin, l, jsmin, jsmax, ll, js, lambda_ket
  INTEGER :: lp1, lp2, lam1, lam2, lambda_bra, iq, iph, lp, is1, is2
  REAL(DP), INTENT(IN) :: pa, pb, pc, pd
  REAL(DP), INTENT(INOUT) :: v_pot
  REAL(DP) :: p_cd, p_ab, q_ab, q_cd, q_rel, k_rel, qmin, qmax, kinetic_ket
  REAL(DP) :: sum6j, sum9j, lsjj_coeff_bra, lsjj_coeff_ket
  REAL(DP) :: int_sum
  REAL(DP) :: vb_bra, vb_ket, kinetic_bra, pos, p_com, vv
  LOGICAL triag
  v_pot = 0.0D0; kinetic_ket =  0.5*(pc*pc+pd*pd)
  kinetic_bra =  0.5*(pa*pa+pb*pb)
  ! Tests of momenta relation needed for the computation of vector brackets
  ! Max values of two-body momenta
  p_cd = pc+pd; p_ab = pa+pb
  ! Min value of two-body momenta
  q_cd=ABS(pc-pd); q_ab=ABS(pa-pb)

  qmin = max(q_cd,q_ab)
  qmax = min( p_cd,p_ab)

  
  q_x(:) = 0.5d0*( qmax+qmin ) + 0.5d0*( qmax-qmin )*x_mesh(:)
  q_w(:) = 0.5d0* ( qmax-qmin ) * wx_mesh(:)
  !Under gir ikke lenger NaN
 ! q_x(:) = (0.5d0*(qmax-qmin)*(t_mesh(:)+2.d0*qmax/(qmax-qmin) -1.d0))
 ! q_w(:) = (0.5d0*(qmax-qmin)*wt_mesh(:))* dsqrt( 1.d0 - t_mesh(:)*t_mesh(:) )
  ! write(*,*) 'hei, en test for debugging' !HER ER DET EN SEGMENT FEIL
  !q_rel = phase* dcmplx(0.5d0*(qmax-qmin)*(t(iq)+2.d0*qmax/(qmax-qmin) -1.d0))
  !dq_rel = phase*dcmplx(0.5d0*(qmax-qmin)*w_t(iq))
  !cheb_weight = dsqrt( 1.d0 - t(iq)*t(iq) )

  !CALL gauss_legendre(qmin,qmax,q_x,q_w,qint_points)
  ! Start performing the transformation for CoM and relative coordinates
  ! to lab frame, ket side first
  is1 = 1; is2=1
  DO spin= 0, 1  ! S = 0 or 1 for nucleons only
     DO js = jmin, jmax

        DO l = abs(js-spin), js+spin
           !DO l=0, lmax  ! relative orbital momentum, (jmax changed to lmax)
           !jsmin=ABS(l-spin); jsmax=l+spin

           IF ( (ABS(isospin_tz) == 1).AND.(MOD(l+spin+ABS(isospin_tz),2)==0)) CYCLE
           !   DO js=jsmin,jsmax    ! Relative angular momentum
           !      IF(js > jmax) CYCLE
           llmin=ABS(j_lab-js); llmax=j_lab+js
           DO ll=llmin,llmax ! CoM orbital momentum
              ! parity test
              IF(MOD(l+ll+lc+ld,2) /= 0) CYCLE    ! parity check ket side
              lambda_min=ABS(lc-ld);  lambda_max=lc+ld
              DO lambda_ket = lambda_min,lambda_max
                 IF(triag(lambda_ket,l,ll)) CYCLE
                 IF(triag(lambda_ket,j_lab,spin)) CYCLE
                 lsjj_coeff_ket = sum6j(ll,l,lambda_ket,spin,j_lab,js)* &
                      sum9j(jc,jd,lc,ld,is1,is2,j_lab,lambda_ket,spin)
                 IF (ABS(lsjj_coeff_ket) < epsilon) CYCLE
                 ! Transformation for CoM and relative coordinates
                 ! to lab frame, bra side now
                 lp1=ABS(js-spin); lp2=js+spin
                 lam1=ABS(j_lab-spin); lam2=j_lab+spin
                 !  write(*,*) lsjj_coeff_ket
                 ! The integral over momenta is the last loop
                 DO iq=1, qint_points

                    p_com = q_x(iq)
                    pos = kinetic_ket-0.25D0*p_com*p_com
                    IF(pos <=  0.D0) CYCLE
                    q_rel = SQrt(pos)
                    if( q_rel > abs(ra(vkk_dim) ) ) CYCLE
                    pos  = kinetic_bra - 0.25D0*p_com*p_com
                    IF ( pos < 0.D0 ) CYCLE; k_rel = SQRT(pos)
                    if( k_rel > abs(ra(vkk_dim) ) ) CYCLE
                    !  write(*,*) kinetic_bra, p_com, pc+pd, DABS(pc-pd)
                    vb_ket = vector_trcoefficients(q_rel,p_com,pc,pd,l,ll,&
                         lambda_ket,lc,ld,isospin_tz)
                    !  write(*,*) 'vb_ket', vb_ket
                    ! TIL HIT ER DET BRA, INGEN UENDELIGHETER I vb_ket
                    IF (ABS(vb_ket) < epsilon) CYCLE !ok
                    int_sum = 0.0_dp
                    DO lp=lp1,lp2
                       IF ((ABS(isospin_tz) == 1).AND.(MOD(lp+spin+ABS(isospin_tz),2)==0)) CYCLE
                       IF(lp > lmax) CYCLE
                       IF(MOD(lp+ll+la+lb,2) /= 0) CYCLE         !  parity check bra side
                       !  get the interaction in the CoM coordinates
                       IF ( ( lp == l) .AND. ( lp == js)) THEN
                          ncoup = 1
                       ELSEIF ( ( lp /= js) .OR. ( l /= js)) THEN
                          ncoup = 2
                       ENDIF
                       CALL get_channel(lp,l,js,spin,isospin_tz,ncoup,channel)
                       IF ( channel == 0) CYCLE


                       !vv_cmplx = 0.0
                       !if ( l > 0 ) cycle
                       !if ( lp /= l ) cycle
                       !                       if ( channel /= 1 ) cycle
                       !                       if ( js > 3 ) cycle

                       !call malfliet_tjon( vv_cmplx, l, phase*k_rel, phase*q_rel )
                       !call gaussian( vv_cmplx, l, js, phase*k_rel, phase*q_rel )

                       !                       if ( l == 0 .and. lp == l ) then
                       !                      call yama( vv_cmplx, k_rel, q_rel )
                       !                       end if



                       CALL g_interpolate(channel, ncoup, lp,l,js,k_rel,q_rel,vv)
                       !CALL v_bare(lp,k_rel,l,q_rel,spin,js,isospin_tz,vv_cmplx)
                       !write(*,*) 'vv_cmplx', vv_cmplx
                       !infinity i vv_cmplx

                       ! don't need, because the vector coeff has taken it into account                    
                       !                       !for complex phase
                       !                      vv=vv*iph(ABS(l-lp)/2)
                       vv=vv*iph(ABS(la+lb-lc-ld)/2)                       

                       if( isospin_tz==0) then
                          vv = vv*iph(ABS(l-lp) )
                       else
                          vv=vv*2.0_dp
                       endif
                       
                       
                       DO lambda_bra=lam1,lam2
                          IF(triag(lambda_bra,la,lb)) CYCLE           ! triangular rels.
                          IF(triag(lambda_bra,lp,ll)) CYCLE
                          lsjj_coeff_bra = sum9j(ja,jb,la,lb,is1,is2,j_lab,lambda_bra,spin)* &
                               sum6j(ll,lp,lambda_bra,spin,j_lab,js)
                          !  write(*,*) 'lsjj',lsjj_coeff_bra ! Ingen uendeligheter
                          IF (ABS(lsjj_coeff_bra) < epsilon) CYCLE
                          vb_bra = vector_trcoefficients &
                               (k_rel,p_com,pa,pb,lp,ll,lambda_bra,la,lb,isospin_tz)
                          !            write(*,*) 'vb_bra', real(vb_bra) ! ingen uendeligheter her
                          IF (ABS(vb_bra) < epsilon) CYCLE
                          !   write(*,*) 'HEI DETTE ER EN TEST'
                          int_sum = int_sum+lsjj_coeff_bra*vv*vb_bra
                          !     write(*,*)  real(vv), real(vb_bra)
                       ENDDO
                    ENDDO
                    !     if ( abs(int_sum) < epsilon ) cycle
                    !  write(*,*)'int_sum', int_sum
                    v_pot=v_pot + int_sum*vb_ket*lsjj_coeff_ket*q_w(iq)*0.25D0
                 ENDDO
              ENDDO

           ENDDO
        ENDDO
     ENDDO
  end DO
  v_pot=v_pot*hbarc**3
END SUBROUTINE v_k

!
REAL(DP) FUNCTION sum9j(j1,j2,l1,l2,is1,is2,jt,lam,is)
  USE constants
  USE ang_mom_functions
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: j1,j2,l1,l2,is1,is2,jt,lam,is
  REAL(DP) :: s1, s2, sls

  s1=SQRT((1.*j1+1.)*(1.*j2+1.)*(2.*is+1.))
  sls=snj(2*l1,2*l2,2*lam,is1,is2,2*is,j1,j2,2*jt)
  s2=2.*lam+1.
  sum9j=sls*s1*s2

END FUNCTION sum9j

REAL(DP) FUNCTION sum6j(lc,l,lam,is,jt,js)
  USE constants
  USE ang_mom_functions
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: lc,l,lam,is,jt,js
  REAL(DP) :: s3, sign, sx

  s3=SQRT(2.*js+1.)
  sign = 1.- 2.*MOD(lam+js+lc+is,2)
  sx=sjs(2*lc,2*l,2*lam,2*is,2*jt,2*js)
  sum6j=sign*s3*sx

END FUNCTION sum6j

SUBROUTINE g_interpolate(channel,ncoup,la,lb,jang,k_rel,q_rel,vv)
  USE wave_functions
  USE constants
  USE relcm_gmatrix
  IMPLICIT NONE
  REAL(DP), INTENT(IN) :: k_rel, q_rel
  REAL(DP), INTENT(INOUT) :: vv
  INTEGER, INTENT(IN) :: ncoup, channel, la, lb, jang
  INTEGER :: lim1, lim2, lim3, lim4, dim

  ! write(*,*) 'n_k, n_k1, n_k2' , n_k, n_k1, n_k2
  vv = 0.0_dp;
  IF( (k_rel < k_cutoff) .AND. (q_rel < k_cutoff) ) THEN
     IF ( ncoup == 1) THEN
        lim1=1; lim2=vkk_dim ; lim3=1 ; lim4=vkk_dim
     ELSEIF ( ncoup == 2 ) THEN
        IF ( (la == lb).AND. ( jang > la) ) THEN
           lim1=1; lim2=vkk_dim ; lim3=1 ; lim4=vkk_dim
        ELSEIF ( (la == lb).AND. ( jang < la) ) THEN
           lim1=1+vkk_dim; lim2=vkk_dim+vkk_dim ; lim3=1+vkk_dim ; lim4=vkk_dim+vkk_dim
        ELSEIF ( la >  lb ) THEN
           lim1=1+vkk_dim; lim2=vkk_dim+vkk_dim ; lim3=1 ; lim4=vkk_dim
        ELSEIF ( la <  lb ) THEN
           lim1=1; lim2=vkk_dim ; lim3=1+vkk_dim ; lim4=vkk_dim+vkk_dim
        ENDIF
     ENDIF

     CALL lagrange_2dim(k_rel,q_rel,REAL(vfree(lim1:lim2,lim3:lim4,channel)), &
          vv,vkk_dim,abs(ra))
  ENDIF
END SUBROUTINE g_interpolate



SUBROUTINE  get_channel(lp,l,js,spin,isospin_tz,ncoup,channel)
  USE constants
  USE partial_waves
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: isospin_tz, spin, js, lp, l, ncoup
  INTEGER, INTENT(INOUT) :: channel
  INTEGER :: loop, la, lb

  channel = 0
  IF ( ncoup == 1) THEN
     la = js; lb = js
  ENDIF
  IF ( ncoup == 2) THEN
     IF ( js > 0) THEN
        la = js-spin; lb = js+spin
     ELSEIF( js == 0 ) THEN
        la = js+spin; lb = js+spin
     ENDIF
  ENDIF
  DO loop=1, no_channels
     IF ((orb_lrel_max(loop) == lb) .AND. (orb_lrel_min(loop) == la) .AND. &
          ( spin_rel(loop) == spin) .AND. &
          (jang_rel(loop) == js) .AND. (iso(loop) == isospin_tz ) ) THEN
        channel = loop
        EXIT
     ENDIF
  ENDDO

END SUBROUTINE  get_channel


!     Obtain the bare interaction in a harmonic oscillator 
!     basis plus the kinetic energy and the Coulomb part. It contains
!     also the CoM correction to the relative coordinates. The latter depends
!     on the mass number of the nucleus
! 
SUBROUTINE vkk_channel(i,ncoup)
  USE wave_functions
  USE relcm_gmatrix
  USE partial_waves
  USE configurations
  USE single_particle_orbits
  USE constants
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: ncoup,i  ! loop variable over NN channels
  INTEGER :: n,np, bra, ket, k1, k2, a, c
  INTEGER :: la, lb, jang, lim1, lim2, lim3, lim4
  REAL(KIND = 8) :: e_coulomb, vsum
  COMPLEX(DPC), ALLOCATABLE, DIMENSION(:,:) ::  vzz(:,:),heff(:,:)
  REAL(KIND = 8), ALLOCATABLE :: vkk(:,:)
  ALLOCATE(vkk (ncoup*n_rel,ncoup*n_rel))
  ALLOCATE(vzz (ncoup*n_rel,ncoup*n_rel))
  vzz = 0.0D0
  jang=jang_rel(i)
  !     setup the NN interaction
  vkk = 0.0D0
  CALL nocorepotential_interface(ncoup,vkk,i)
  vzz = vkk
  WRITE(6,'(12H Channel Nr:,I3,7H l_min:,I3,7H l_max:,I3,3H J:,I3,3H S:,I3,4H Tz:,I3)') &
       i, orb_lrel_min(i), orb_lrel_max(i), jang_rel(i), spin_rel(i), iso(i)

  ALLOCATE(heff (ncoup*vkk_dim,ncoup*vkk_dim))
  heff =0.0D0
  IF ((type_of_pot == 'OPEP').OR.(type_of_pot == 'Tensorinteraction')   &
       .OR.(type_of_pot == 'LSinteraction') ) THEN
     heff = vzz
  ELSEIF (type_of_renormv =='vbare') THEN
     heff = vzz
  ELSE
     ! get similarity transformed interaction
     CALL vlowk_mtx(ncoup,vzz,heff,i)
  ENDIF
  vfree(1:ncoup*vkk_dim,1:ncoup*vkk_dim,i) = heff(1:ncoup*vkk_dim,1:ncoup*vkk_dim)
  DEALLOCATE(vkk); DEALLOCATE(vzz,heff)
END SUBROUTINE vkk_channel

!
!                 Set up h.o. wf for rel system
!
SUBROUTINE ho_wfk
  USE constants
  USE wave_functions
  IMPLICIT NONE
  INTEGER :: n, l, i, j
  REAL(KIND = 8)  :: cx(0:200), factor, z_rel, xp, ph, oscl_r, sum_rel, contrib
  
  DO n=0,nmax
     ph=(-1.D0)**n
     DO l=0,lmax
        sum_rel=0.
        factor = 0.5D0*((n+l+2)*LOG(2.D0)+fac(n)-dfac(2*n+2*l+1)-0.5D0*LOG(pi))
        factor = EXP(factor)
        DO i=1,n_lab
           z_rel= k_mesh(i)*k_mesh(i)*oscl*oscl
           CALL laguerre_general( n, l+0.5D0, z_rel, cx )
           xp = EXP(-z_rel*0.5)*((k_mesh(i)*oscl)**l)*cx(n)
           rnlr(i,l,n) = xp*(w_mesh(i)*k_mesh(i))*ph*factor*(oscl**(1.5D0))     ! rel wf
           contrib = xp*ph*factor*(oscl**(1.5D0)) 
           sum_rel=sum_rel+ w_mesh(i)*(contrib*k_mesh(i))**2
        ENDDO
!        WRITE(*,*) 'Norm rel ho wf n,l : ', n, l, sum_rel
     ENDDO
  ENDDO

END SUBROUTINE ho_wfk
